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16 Abstract 

A method for calculating transonic flow over steady and oscillating airfoils has 
been developed by Isogai. It solves the full potential equation with a semi-implicit, 
time-marching, finite difference technique. Steady flow solutions are obtained from 
time asymptotic solutions for a steady airfoil. Corresponding oscillatory solutions 
are obtained by initiating an oscillation and marching in time for several cycles 
until a converged periodic solution is achieved. In this paper the method is 
described in general terms and results for the case of an airfoil with an oscillating 
flap are presented for Mach numbers 0.500 and 0.875. Although satisfactory results 
are obtained for some reduced frequencies, it is found that the numerical technique 
generates spurious oscillations in the indicial response functions and in the 
variation of the aerodynamic coefficients with reduced frequency. These oscillations 
are examined with a dynamic data reduction method to evaluate their effects and 
trends with reduced frequency and Mach number. Further development of the numerical 
method is needed to eliminate these oscillations. 
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INTRODUCTION 


The flutter critical portion of the aircraft flight envelope generally 
occurs at transonic speeds. This critical condition results from both the 
high dynamic pressures of operation and the dip in flutter speed or "bucket" 
that occurs at transonic Mach numbers. The dip in flutter speed is influenced 
by airfoil thickness and shape and cannot be satisfactorily treated by state- 
of-the-art aerodynamic analyses. Thus an important current topic in aero- 
dynamic research is the development of methods for the calculation of unsteady 
aerodynamics for use in transonic flutter analysis. Many of the current 
efforts have built on the recent success of steady flow numerical finite 
difference solution procedures and, to date, have primarily been applied to 
two-dimensional airfoils as a means of evaluating and refining the analyses 
and algorithms involved. One such method of solving the full potential 
equation was presented in reference 1. It has been extended by Isogai 
(unpublished) and has been applied to calculate several examples of unsteady 
transonic flow (references 2-3). In general, comparisons with other methods 
were favorable (references 1-2). One exception was that strong shocks were 
located further aft on the airfoil than corresponding Euler equation solutions 
(ref. 1). This trend might be anticipated since the full potential equation 
overpredicts the strength of shock waves (ref. 4). Similar comparisons with 
experimental data of references 5 and 6 for the NACA 64A006 airfoil with an 
oscillating flap were also in reasonable agreement (references 2-3) considering 
the fact that the theory did not allow for viscous or wind tunnel effects. 

To define the variation of aerodynamic forces with reduced frequency for 
a transonic flow, the forces were calculated using Isogai' s method for a 
moderate set of k-values and for a Mach number of 0.875 (ref. 3). Results 
of these calculations showed oscillations in the variation of the forces with 
k which were particularly pronounced for reduced frequencies less than 0.10. 

In this paper a more extensive set of calculations are presented for a Mach 
number of 0.500 along with the previous results for a Mach number of 0.875. 

The oscillatory behavior of the coefficients with k is shown to exist also 
at the lower Mach number. The oscillations are believed to be spurious and 
generated by the numerical method. Corresponding indicial response functions 
for flap displacement are also presented.. They are analyzed by the dynamic 
data reduction method of reference 7 in an effort to define further and to 
evaluate the effects of the oscillations and their trends with Mach number. 
Other calculations are also discussed in which the finite difference grid 
and the time step were varied in order to define more precisely the source of 
the oscillations. 


NOMENCLATURE 


a (speed of sound) /V 

a^ (freestream speed of sound) /V 
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imaginary or out-of-phase part 
reduced frequency, toc/2V 


reduced frequency of critical or "modal" oscillation 


Mach number 
pressure (dimensional) 

2 

freestream dynamic pressure Jjp V 

real or in-phase part 
nondimens ional time, t/(c/2V ) 

OO 

dimensional time, s 

nondimensional time step 

(velocity component in x-direction)/^ 

(velocity component in y-direction)/V Qo 

free-stream total velocity (dimensional) 


coordinate distances 
angle of attack 


positive 


nose up. 
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y ratio of specific heats 

flap deflection, positive trailing edge down 

0 pitch angle, positive nose up 

f ree-stream density (dimensional) 

cj) perturbation potential 

a) frequency of oscillation, rad/s 

GENERAL DESCRIPTION OF THE METHOD 


The potential equation for two-dimensional time-dependent flow is 

(a 2 - U 2 )4> - 2UV(J) + (a 2 - V 2 )<{> 

xx xy yy 

- 2U<j> - 2V<j> „ - <j> = o 

xt yt tt 


where 


(J) = perturbation potential 

U = cos a + d) 

x 

V = sin a + (j) 

y 

a 2 = a 2 - 0.5(y - 1) (24> t + U 2 + V 2 -1) 

This equation is nonlinear because a, U, and V are functions of cf> , and 
numerical finite difference techniques are generally used to obtain solutions. 
Use of the full potential equation yields a method that is intermediate in 
completeness and computational effort between methods that use the Euler 
equations (ref. 8) and those that use the small disturbance equation (references 
9 and 10). One advantage of using the potential equation as compared with the 
Euler equations is that for the potential equation only the single variable 
<(> has to be stored, whereas for the Euler equations p, p, U, and V must 
be stored. 

The semi-implicit finite difference technique (ref. 1) is used which 
alleviates some of the limits on time step required for stability. An unpub- 
lished addition to the algorithm of reference 1 has also been made by Isogai 
which increases the permissible time step. It essentially is an additional 
implicit pass through the flow field. The time step is still restricted, 
however, which implies that the computer time increases as reduced frequency 
k decreases. The method also uses a rotated difference scheme to maintain 
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numerical stability. 

A stretched rectangular Cartesian grid system similar to that of reference 
11 is used to map the infinite physical space to a finite computational system. 
For all calculations presented herein, a 57 x 57 grid system was used, although 
some calculations using other grid systems are discussed. The upper half of 
the grid is shown in figure 1 in physical space. One finite y-grid line and 
the lines at infinity are not shown. The mapping clusters points near the 
leading and trailing edges and near the airfoil surface, but no special consid- 
eration is given to shocks or hinge lines. For the grid system used (fig. 1), 
there are 29 points on each airfoil surface with the first and last points 
on the airfoil at 0.01c and 0.99c, respectively. In comparison with conven- 
tional steady flow transonic calculations, this grid system might be considered 
a medium grid rather than a coarse or fine grid system. 

The airfoil motion boundary condition is applied at the mean airfoil 
position. This assumption may restrict the valid range of amplitude of 
motion but simplifies the computer program because the airfoil and computation 
grid remain fixed in time. In addition, the finite difference method uses 
a quasi-conservative shock-capturing difference scheme to treat moving shock 
waves (ref. 1). This treatment of shock waves ensures that the correct shock 
jump relations are maintained. It also simplifies the computer program since 
the moving shocks are treated automatically. They are, however, smeared over 
a few mesh spaces. 

The current version of the computer program is dimensioned for a grid of 
up to 61 x 61 points. It requires 37 000 (114 Kg) locations of central 
computer memory. On a CYBER 175, operating under the NOS operating system 
with the FTN compiler, the program takes about % seconds of CPU time for each 
time step. About 3000 time steps are required to converge for a steady flow 
case and 1000 to 5000 steps for a typical oscillatory case. 

RESULTS AND DISCUSSION 
Calculations for M = 0.500 

The calculated pressure coefficient for the 64A006 airfoil in steady 
flow at zero angle of attack is compared with experimental data from reference 
5 in figure 2. The theory is generally in good agreement but there is some 
small deviation from the experimental data over the aft portion of the airfoil. 

Calculations for harmonic motion .- The variation of unsteady lift with 
reduced frequency k for a quarter-chord flap oscillating with an amplitude 
of 1.5° is shown in figure 3 as calculated from linear theory and from the 
present finite difference method using the nonlinear full potential equation. 
Harmonic flap oscillation was assumed, and a time step At = .05 was used. 

For this subsonic case there should be good agreement with linear theory with 
the possible exception of near k = 0. Although the underlying trend of the 
nonlinear theory is similar to that of linear theory, there are several severe 
excursions in the results that are apparently spurious. It is believed that 
these excursions are introduced by the numerical technique. Examination of 
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the output of the program indicates that, at all points in the flow field, 
there are oscillations which are in phase. Thus the oscillations appear 
more nearly like standing waves than traveling waves. At reduced frequencies 
near the peaks of the oscillations, the unsteady load distribution is shifted 
in phase near the leading edge such that the primary influence is on c^. 

Here c, is shown as it is a more sensitive indicator than c or c, . 

1 m h 

The two theories also give a somewhat different behavior near k = 0 
where thickness may be more important. However, for the finite difference 
calculations at k = 0, the circulation has been neglected for the far field 
boundary in all of the calculations of this report. In steady flow, the 
far field boundary should include the mean circulation. If the circulation 
is included in the calculations for figure 3, the coefficients increase by 
about 10%. (Such a result may suggest that the finite difference grid may 
not be adequate.) Although the far field circulation can be readily included 
in calculations for k = 0, no similar boundary condition is available for 
time dependent calculations, and it is consistent to neglect it for both 
steady and unsteady calculations. Note also that the trend of the two 
methods near k = 0 may indicate a somewhat different behavior near k = 0, 
but further definition is required to obtain a conclusive result. 

Indicial calculations .- A large number of points were required to define 
the trends shown in figure 3. The 52 points shown and the steady flow 
calculation required a significant expenditure of computer resources, 
particularly for the points at low reduced frequencies. The same finite 
difference program was also used to generate the response for indicial flap 
deflection shown in figure 4. The response has the expected overall shape 
of a jump at t = 0 plus a real exponential growth to a steady state value. 

In addition there are several oscillations of different character that appear 
in the response. Both this curve and the appearance of the results of 
figure 3 suggest that the response is somewhat analogous to the response of a 
multi-mode mechanical system and that such an analogy might be useful in 
defining the characteristics of the results. The method of reference 7 has 
been applied to fit the response of figure 4 with a series of complex expo- 
nential functions in a least square sense such that the amplitude, phase, 
frequency, and damping of each "mode" or term in the series is calculated. 

The fit using 8 terms is shown to be essentially coincident with the transient 
response in figure 4. The components of the transient lift determined by 
the 8-term fit are shown in figure 5. The offset or steady state value and 
the first mode are shown in figure 5a and the oscillatory terms are shown in 
figure 5b. For a steady state calculation, the rate of approach to convergence 
would be controlled by mode 2, and there should be no influence on the 
results after it dies out. For the time-dependent problem, however, these 
"modes" significantly distort the unsteady forces as shown in figure 3. A 
further indication of these effects is given in figure 6. The harmonic 
calculations (from fig. 3) are compared with the lift coefficient calculated 
from the analytic Fourier transform of the 8-mode fit. Good agreement is 
obtained with some differences in the higher range of k-values . Calculation 
of the forces for harmonic motion from the Fourier transform of the indicial 
response is valid only for linear systems. 
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The response for indicial flap displacement should have a delta function 
at t = 0+ that results from the discontinuous velocity at t = 0. This is 
in contrast to the usual indicial functions for subsonic compressible flow 
that are traditionally given in terms of an indicial velocity rather than 
displacement (ref. 12). The delta function for the flap displacement contri- 
butes to the unsteady lift at all k-values by 

Ac l = 1 4M (1) 

Equation 1 has been derived in a manner similar to the discussion of 
reference 12. The results of the transient analysis shown in figure 6 include 
the addition of this term to the transform of the 8 term fit. The finite 
difference scheme takes a finite time step at t = 0 and thus cannot resolve 
the details of the transient sufficiently to include the delta function. Note 
that equation (1) gives a large contribution to the imaginary part of c 1 (1.0 
at k = 2.0 and M = 0.5). 

The use of the indicial function and the data reduction technique of 
reference 7 give essentially the same information as that contained in figure 
3 but at a cost comparable to the cost of 2 or 3 low frequency points of 
figure 3, a substantial savings. It appears that potential problems with 
difference techniques can be analyzed economically in this manner. The 
additional variable of time makes a global verification of a program of this 
type a significantly larger task than that for a similar steady flow program. 

Relationship to numerical stability .- The question arises as to the 
relationship of the oscillations to the time step limited stability criterion. 
A transient for a marginally unstable case is shown in figure 7 which is for 
a At = 0.09825, about twice the At = 0.05 used previously. The first 
portion of the curve is qualitatively very similar to that of figure 4 
(allowing for the difference in time scale). The latter part of the curve, 
however, (fig. 7) shows a very rapid divergence with some small high frequency 
oscillations superimposed. Thus it appears that the modes previously dis- 
cussed are not those involved in the numerical instability, but rather a 
static "divergence" is involved. In contrast, the difference scheme of 
reference 1 goes unstable in an oscillatory fashion at 2/3 of the Nyquist 
frequency, k = 2/3 (^-) . The additional smoothing pass of the present method 

apparently suppresses the high frequency instability and thus increases the 
allowable time step. 

Effects of Various Parameters on Oscillations 

To define further some aspects of the oscillations, parameters such as 
time step, finite difference grid, and Mach number have been explored. As 
discussed previously, much of the transient response at M = 0.5 was 
essentially the same for At = 0.05 and for At = 0.09825. This situation 
was also the case for smaller time steps except very near t = 0, and the 
results appear to be remarkably insensitive to time step. 
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Calculations have been made with the number of points used in the finite 
difference grid reduced by half in the x-direction, half in the y-direction, 
and by half in both directions. Reducing the number of points in the x- 
direction had only a modest influence on the frequencies of the oscillations. 
However, with only half the points in the y-direction, the frequencies of 
the oscillations were approximately twice the frequencies obtained with the 
full grid. Reducing the number of points in both directions also doubled 
the frequencies of oscillations. It appears therefore that the present 
results are a strong function of the y-grid, and further investigation is 
needed. 

/■ 

The transonic perturbation method of reference 9 has had considerable 
difficulty with frequency limitations that have prevented meaningful solutions 
using conventional relaxation methods except for low values of k. Recently, 
it has been found that by solving the large system of simultaneous linear 
equations with an out-of-core solver, good results can be obtained. Frequency 
limitations can still be encountered, but they are very sensitive to the 
y-grid arrangement and can be alleviated by the proper choice of finite 
difference grid. This finding is consistent with the present results for a 
time dependent method. 

Calculated results for a 1° pitching oscillation are shown in figure 8 
for a range of reduced frequency up to k = 0.5. Similar deviations occur 
at the same reduced frequency as for the oscillating flap. 

The frequencies of oscillations k^ have been determined for several 

Mach numbers by using the indicial response, inspecting the Fourier transform 
for peaks and by fitting with complex exponentials in a manner similar to the 
calculations M = 0.5. The resulting frequencies are shown in figure 9. 

For the lower Mach numbers, the frequencies follow an orderly pattern and 
decrease with Mach number. The damping of the oscillations tends to increase 
with Mach number, particularly for the higher frequencies, so that there are 
relatively fewer identifiable frequencies above M = 0.8. Furthermore, in 
the transonic range for M = 0.85 or higher, shock waves tend to invalidate 
the linear assumptions. For M = 0.8 to 0.9, the method may be limited to 
values of k larger than the values of k c shown as symbols in figure 9. 

Calculations for M = 0.875 

Steady flow calculations .- The calculated steady pressure coefficient 
for the 64A006 airfoil is compared with experimental values (ref. 6) in 
figure 10. The calculated shock wave is stronger and further aft than the 
experimental one and is fairly near the hinge line at .75c. 

The corresponding flap hinge moment for static deflections, shown in 
figure 11, indicates a strong nonlinearity. As flap deflection is increased, 
the slope of the curve is initially positive (unstable) , rapidly changes to a 
large stable value, and then levels out. This nonlinear behavior is brought 
about by the movement of the shock over the upper surface of the flap as the 
flap deflection is increased and is illustrated by the pressure distributions 
of figure 12(a)-12(c). Initially, as a result of the relative movement of 
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the upper and lower surface shock waves, there is a small down load over the 
flap (fig. 12a). At an intermediate deflection level, there is a large 
effect on the loading of the flap and airfoil as the upper surface shock moves 
aft over the airfoil and the lower one moves forward (fig. 12b). However, 
once the upper surface shock reaches the trailing edge, there is decreased 
effectiveness as only the loading over the flap is changed. Although this 
trend may not be accurately predicted by the present method (since shock 
wave strength is overpredicted) such a phenomenon would be anticipated to 
occur for conditions that lead to movement of a shock over the control 
surface. In practice, however, these trends might also be altered signif- 
icantly by boundary layer effects. This case illustrates a phenomenon 
that should be interesting to investigate more fully, particularly consider- 
ing its aeroelastic implications. 

Unsteady calculations .- The transient lift for an indicial flap displace- 
ment of 1.5° corresponding to this nonlinear static behavior is shown in 
figure 13. There is a considerable change in character from that at 
M = 0.5 (fig. 4). The initial portion is nearly linear and then a spike or 
corner occurs, apparently as a result of shock motions. The latter portion 
of the response again involves an oscillatory behavior. The exponential 
series fits the response well except near the corner or spike. 

Oscillatory results calculated from harmonic motion are compared with 
linear theory in figure 14. The nonlinear results approach the linear 
results at the higher k-values, but there are again large oscillations of the 
unsteady lift with reduced frequency in the results of the finite difference 
calculations. Although the results would not be expected to agree with linear 
theory at low reduced frequencies, these oscillations are not considered to 
be physical but introduced by the numerical method. In particular, the result 
at k = 0.059 is near one of the frequencies extracted from the transient 
response. It might be noted that these results at M = 0.875 smooth out much 
more rapidly with k (fig. 14) than do the results of M = 0.5 (fig. 3). 

The unsteady coefficients have also been calculated from the Fourier 
transform of the computed fit to the indicial response and are compared 
with the calculations for harmonic oscillations (1.08° flap amplitude) 
in figure 15. Although the two results are somewhat similar in character, 
they do not agree well at low reduced frequencies. The transform process 
depends on linearity and the static forces are quite nonlinear in this case 
(fig. 11). The indicial response and its tranform however would be a good 
indicator of possible problems with the difference method. 

CONCLUDING REMARKS 

Extensive calculations with a semi- implicit method for solving the two- 
dimensional unsteady full potential equation have indicated that the numerical 
technique as presently formulated generates spurious oscillations in the 
indicial response functions and in the corresponding aerodynamic coefficients 
versus reduced frequency. These oscillations appear to be related to the 
spatial finite difference grid, and further development is needed to improve 
the accuracy of the method. 
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An extensive computational effort is required to verify a time dependent 
program such as the one discussed herein, but this effort can be reduced 
through the use of indicial methods. In the present investigation, the use 
of conventional dynamic data analysis methods of Fourier transforms and of 
time-plane fitting with complex exponentials was helpful in determining the 
existence and trends of the oscillations. 
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Figure 1.-- Grid system used in calculations, 
57 x 57. 



Figure 2.- Steady pressure distribution for the 
NACA 64A006 airfoil at M = 0.5 and a = 0. 
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HARMONIC OSCILLATION 
TRANSIENT ANALYSIS 



Figure 6.- Comparison of harmonic and transient 
analyses for lift due to flap motion at M = 0.5 

and a = 0. 



Figure 7.- Transient lift for 1.5° step in flap 
deflection at M = 0.5 and for At = 0.09825. 
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Figure 12.- Steady pressure distribution for the 
NACA 64A006 airfoil at M = 0.875 and a = 0. 

(a) S f = 0.15° 


Figure 12.- Continued, 
(b) 6 f = 0.60° 





Figure 12.- Concluded. Figure 13.- Transient lift for a 1.5° step in flap 

deflection at M = 0.875 and for At « 0.11. 

(c) 6 f = 1 .08° 



k 

Figure 14.- Unsteady lift due to flap motion at 
M = 0.875 and a = 0. 



Figure 15.- Comparison of harmonic and transient 
analyses for lift due to flap motion at M = 0.875 
and a = 0. 
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